#include <TROOT.h>
void Z_CHIRALSEA(char *momentum_order, int beta100, int l, int t, int nclust, char *dir, int imom, char *pmprint, char *Oa2print, char *rgiprint , int iprop)
{

        L=l;
        T=t;
        Nclust=nclust;

        double ZAPREV=0.757;
        int Nmass;
        int set_p[4];
        char setpprint[16];
 


	  // ensemble-mu-mpcac input file

	  char mu_dat[1024];
	  sprintf(mu_dat,"../INPUT/masses_%d.dat",beta100);
	  FILE *mudat=fopen(mu_dat,"r");
	  fscanf(mudat,"%i",&Nmass);

	
	  char letter[8];
	  double *mu=(double*)malloc(sizeof(double)*(Nmass));
          double *mPCAC=(double*)malloc(sizeof(double)*(Nmass));
	  int *e=(int*)malloc(sizeof(int)*(Nmass));
	  fscanf(mudat,"%s",letter);
          for (int im=0; im<Nmass; im++)
                fscanf(mudat,"%d",&e[im]);
	  for (int im=0; im<Nmass; im++)
		fscanf(mudat,"%lf",&mu[im]);
	  for (int im=0; im<Nmass; im++)
		fscanf(mudat,"%lf",&mPCAC[im]);
	  fscanf(mudat,"%lf",&ZAPREV);
	


          //// Momenta set
    whichProp(iprop,setpprint);

    char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
    read_momenta_data(momentafile);
    double *nn=(double*)malloc(sizeof(double)*(Nmom));
    compute_nn(momentum_order, nn);
    if (iprop<16)
    {
    whichSetp(iprop,set_p);
    Momentum_Magro(momentum_order, set_p);
    }
    else
    {
    Momentum(momentum_order);
    }
    // equivalent momenta
    Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
    Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);
	printf("Nmom eq %d\n",Nmom_eq);



     //////////////////////////////////
     // IO files
     ///////////////////////////////
     char file[1024];
     char z_output[1024];
     char gammaprint[16];



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

  
	//////////////////////////////////
	// Read Zs
	/////////////////////////////

         double **Zq=malloc(sizeof(double*)*(Nclust+1));
         for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double)*(Nmass));

         double ***ZO=malloc(sizeof(double**)*(5));
         for (int igamma=0; igamma<5 ; igamma++) ZO[igamma]=malloc(sizeof(double*)*(Nclust+1));
         for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust]=malloc(sizeof(double)*(Nmass));

	double *aux_clust=malloc(sizeof(double)*(Nclust+1));


    // Zq

    for (int imass=0; imass<Nmass; imass++)
    {
		if ( strcmp(letter,"NULL")==0)  	sprintf(file,"%s/%d%s/%s/QAV/p%03d/Zq_%s%s.jk",dir,e[imass],pmprint,setpprint,imom,Oa2print,rgiprint);
		else sprintf(file,"%s/%s%d%s/%s/QAV/p%03d/Zq_%s%s.jk",dir,letter,e[imass],pmprint,setpprint,imom,Oa2print,rgiprint);
//		printf("Reading %s\n", file);
                FILE *fin=fopen(file,"r");
                if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",file);
                 }

                 for (int iclust=0; iclust<Nclust+1; iclust++) fscanf(fin,"%lf",&aux_clust[iclust]);
                 fclose(fin);

	 for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imass]=aux_clust[iclust];
    }

    // ZO 
   
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        {
	   for (int imass=0; imass<Nmass; imass++)
	   {

	  	if ( strcmp(letter,"NULL")==0)  sprintf(file,"%s/%d%s/%s/QAV/p%03d/Z%s_%s%s.jk",dir,e[imass],pmprint,setpprint,imom,gammaprint,Oa2print,rgiprint);
		else sprintf(file,"%s/%s%d%s/%s/QAV/p%03d/Z%s_%s%s.jk",dir,letter,e[imass],pmprint,setpprint,imom,gammaprint,Oa2print,rgiprint);
                FILE *fin=fopen(file,"r");
                if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",file);
                 }

                 for (int iclust=0; iclust<Nclust+1; iclust++) fscanf(fin,"%lf",&aux_clust[iclust]);
                 fclose(fin);

         for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust][imass]=aux_clust[iclust];

	 }

       }
    }



      ///////////////////////////////////////
      // CHIRAL FIT FOR EACH MOMENTUM and PM
     ///////////////////////////////////////

      // auxiliar for minuit routines
      double *ZAUX=malloc(sizeof(double)*(Nmass));
      double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmass));
      for (int imass=0; imass<Nmass; imass++) ZAUX_CLUST[imass]=malloc(sizeof(double)*(Nclust+1));

      double *MAUX=malloc(sizeof(double)*(Nmass));

       double **ZCHI=malloc(sizeof(double*)*(6)); // One extra for Zq
       for (int igamma=0; igamma<6 ; igamma++) ZCHI[igamma]=malloc(sizeof(double)*(Nclust+1));

       double **invZCHI=malloc(sizeof(double*)*(6)); // One extra for Zq
       for (int igamma=0; igamma<6 ; igamma++) invZCHI[igamma]=malloc(sizeof(double)*(Nclust+1));

       double **ZLIN=malloc(sizeof(double*)*(6)); // One extra for Zq
       for (int igamma=0; igamma<6 ; igamma++) ZLIN[igamma]=malloc(sizeof(double)*(Nclust+1));


       double **chi2=malloc(sizeof(double*)*(6)); // One extra for Zq
       for (int igamma=0; igamma<6 ; igamma++) chi2[igamma]=malloc(sizeof(double)*(Nclust+1));
	

     // Chiral fit for Zq
    
    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++) 
     {
                MAUX[imass]=pow(mu[imass],2)+pow(ZAPREV*mPCAC[imass],2);
                ZAUX[imass]=Zq[iclust][imass];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=Zq[iclust_loop][imass];
     }
     lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,ZCHI[5][iclust],ZLIN[5][iclust],chi2[5][iclust]); // Zq is saved in the last element (igamma=5)
        sprintf(file,"%s/%s/%s/QAV/p%03d/Zq%s_%s.ext",dir,pmprint,setpprint,imom,rgiprint,Oa2print); 
	write_output_sea_chiral_extrapolation(file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,ZCHI[5],ZLIN[5],chi2[5]);

     }
	
    

     for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) 
     {	
        whichGammaProm(igamma,gammaprint);
	// Polar mass ^2:
       for (int imass=0; imass<Nmass; imass++)
       {
		MAUX[imass]=pow(mu[imass],2)+pow(ZAPREV*mPCAC[imass],2);
                ZAUX[imass]=ZO[igamma][iclust][imass];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=ZO[igamma][iclust_loop][imass];
       } 
	lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,ZCHI[igamma][iclust],ZLIN[igamma][iclust],chi2[igamma][iclust]);
        sprintf(file,"%s/%s/%s/QAV/p%03d/Z%s%s_%s.ext",dir,pmprint,setpprint,imom,gammaprint,rgiprint,Oa2print);
	write_output_sea_chiral_extrapolation(file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,ZCHI[igamma],ZLIN[igamma],chi2[igamma]);

     }


    // write output

   // Zq
    

	printf("Written output in %s/%s/%s/QAV/p%03d/Z*_%s%s.dat/jk\n",dir,pmprint,setpprint,imom,Oa2print,rgiprint);

   for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZCHI[5][iclust];
   sprintf(file,"%s/%s/%s/QAV/p%03d/Zq_%s%s",dir,pmprint,setpprint,imom,Oa2print,rgiprint);
   write_output_fixedmomenta_real(file,aux_clust,SinP2_eq,imom);
  


  // ZO
  

   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
	for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZCHI[igamma][iclust];
        sprintf(file,"%s/%s/%s/QAV/p%03d/Z%s_%s%s",dir,pmprint,setpprint,imom,gammaprint,Oa2print,rgiprint);
        write_output_fixedmomenta_real(file,aux_clust,SinP2_eq,imom);

    }

    delete Zq; delete ZO;
    delete aux_clust;
    delete ZAUX; delete ZAUX_CLUST; delete MAUX;
    delete ZCHI; delete invZCHI; delete ZLIN;
  

}
